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Many-body theory is largely based on self-consistent equations that are constructed in terms of 
the physical quantity of interest itself, for example the density. Therefore, the calculation of impor¬ 
tant properties such as total energies or photoemission spectra requires the solution of non-linear 
equations that have unphysical and physical solutions. In this work we show in which circumstances 
one runs into an unphysical solution, and we indicate how one can overcome this problem. More¬ 
over, we solve the puzzle of when and why the interacting Green’s function does not unambiguously 
determine the underlying system, given in terms of its potential, or non-interacting Green’s func¬ 
tion. Our results are general since they originate from the fundamental structure of the equations. 

The absorption spectrum of lithium fluoride is shown as one illustration, and observations in the 
literature for some widely used models are explained by our approach. Our findings apply to both 
the weak and strong-correlation regimes. For the strong-correlation regime we show that one cannot 
use the expressions that are obtained from standard perturbation theory, and we suggest a different 
approach that is exact in the limit of strong interaction. 

PACS numbers: 71.10.-w, 71.15.-m, 31.15.ee 


In condensed-matter physics, important formalisms for 
predicting and understanding material properties, such 
as density-functional theory (DFT), many-body pertur¬ 
bation theory (MBPT), and dynamical mean-field theory 
(DMFT), avoid the use of the full many-body wave func¬ 
tion and are instead based on simpler quantities, such as 
densities and Green’s functions. The one-body Green’s 
function G, for example, yields expectation values of one- 
body operators, and the total energy. In particular, G 
gives access to photoemission spectra, the most direct 
experimental observation of electronic structure. G is 
well defined as expectation value of particle addition and 
removal to the Al-body ground state or thermal equilib¬ 
rium. However, to use this definition one should know 
the many-body wave function, which is out of reach. A 
framework where the many-body wave function does not 
appear explicitly is provided by many-body perturbation 
theory [l|, where the interacting Green’s function G is 
given as a functional of the non-interacting Green’s func¬ 
tion Go and the bare Coulomb interaction Vc- An im¬ 
portant idea of MBPT is to avoid a possibly ill-behaved 
perturbation expansion of G in terms of Vc and Gq using 
Dyson equations. These are integral equations that de¬ 
scribe the propagation of particles in terms of an effective 
potential or interaction. For the one-body Green’s func¬ 
tion, for example, this effective potential, which is the 


kernel of the Dyson equation, is the self-energy E. The 
power of this approach resides in the fact that even a 
low-order approximation for E yields contributions to all 
orders in Vc ■ Following Luttinger and Ward Q , E is usu¬ 
ally expressed as a functional of G instead of Gq. How¬ 
ever, this makes the Dyson equation non-linear, which 
leads to multiple solutions. [3|, Ij 

This is a very fundamental and general problem. It 
is different from usual convergence problems, which are 
readily detected, for example from the oscillatory behav¬ 
ior of the results; see the supplemental material for an 
examjAe of such a problem in the case of Hartree Fock 
(HF) [^. The appearance of fully converged, but unphys¬ 
ical results, instead, is much more subtle and dangerous, 
and it has important consequences. It is the main topic 
of the present work. 

We show that the presence of multiple solutions has an 
impact that reaches far beyond numerical problems, and 
even points to cases where the currently used strategies 
to derive approximations break down. We develop our 
ideas using a model that represents the structure of the 
exact theory. Galculations for a real system and links 
to recent observations [1,0 demonstrate the potential of 
this approach for analysis and prediction. In particular, 
we answer several very general and important questions: 
i) Is the problem of multiple solutions specific for cer- 


tain cases, or is it a fundamental problem? Does it only 
appear for Green’s functions or also in the framework 
of other well-established and widely used methods, like 
density-functional based approaches? ii) Does it impact 
calculations for real materials?; Hi) How can one detect 
and avoid unphysical solutions?; iv) Does the problem 
depend on the interaction strength, and are there con¬ 
sequences for many-body theory of strongly interacting 
materials? 

To analyze the problem we use the so-called one-point 
model (0PM) [§-^. This model is not system specific 
and that can be solved exactly, such that the physical so¬ 
lution is well defined. It represents important structural 
aspects of the many-body problem, while collapsing all 
arguments of the Green’s functions, self-energy, and the 
interaction to one point, making the equations scalar. In 
Ref. Q, an approximate version of the 0PM was used 
to discuss multiple solutions within the framework of the 
GW approximation [ll| to the self-energy. 

In the present work we use the 0PM without approx¬ 
imations, which simulates the full many-body problem. 
The exact 0PM Green’s function was derived in Ref. 121 
from the one-point equivalent of the equation of motion 
of G, expressed as a functional differential equation [l^ . 
The exact solution reads 


y[yo,u] = 


yo 


1 -b iuyg 


and s[yo,u] =-^uyo, (1) 
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FIG. 1. (Color online) One-point model (OPM): as a 

function of the interaction V. Squares (red): Zq and solution 
of scheme (A); circles (blue): Z^ and solution of scheme (B); 
continuous line (orange): the exact solution Yq. Inset: 
as a function of the interaction V. 


While neither of the two schemes has convergence prob¬ 
lems when iterating, only scheme (I) converges to the 
physical solution, whereas scheme (II) converges to the 
unphysical solution. This happens because the itera¬ 
tion leads to the continued fraction representation of the 
square root Q in Eq. 


where y, yo, and u represent G, Gq, and Vc, respectively. 
The self-energy s is determined from the Dyson equation 
s[yoi u] = Vo^ ~y~^{yow]- Here s is given as a functional 
of the bare interaction u and the noninteracting Green’s 
function yo- Usually, however, one works with the self¬ 
energy given as a functional of the dressed Green’s func¬ 
tion, s[y,u]. Then the Dyson equation reads 

y = yo + yQs[y,u]y. (2) 


This is, in general, a non-linear equation. We first con¬ 
sider the HF self-energy, which in the OPM is s^^[y, u] = 
— \uy. Let us look at the map Gq —>■ G, i.e., the usual 
case, where yo is set by the system, and one searches y. 
The Dyson equation has two solutions, 


y± — _ 

^HF - y 


-l±Vl + 2V 


( 3 ) 


with the rescaled quantities V = y/yo and V = uy^. Here 
is the physical solution, since it connects smoothly 
to To = 1 at U = 0, and is an unphysical solution, 
that diverges for vanishing interaction. Both are shown 
in the inset of Fig. [TJ In real problems Dyson equations 
are solved iteratively. Two possible iteration schemes are: 
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for X = 2V. The sign of the square root is determined by 
the continued fraction in the iterative procedure 141. 

The simple but general structure of the OPM suggests 
that the same picture should emerge for any Dyson-like 
equation. For example, optical properties and screening 
can be calculated by solving the Bethe-Salpeter equation 
(BSE) for the two-particle correlation function, using the 
GW approximation for the self-energy [l^ . The screened 
Coulomb interaction W is calculated from the BSE, and 
is also part of its kernel. Like in the above HF case, this 
makes the problem in principle self-consistent (see e.g. 
Ref. [13). Alternatively, one can use time-dependent 
density-functional theory (TDDFT), that obeys a sim¬ 
ilar Dyson-like equation for the reducible polarizability 
X III- Here we use TDDFT within the so-called boot¬ 
strap approximation proposed in Ref. [l^ . The corre¬ 
sponding Dyson equation for frequency w reads 


X{^) = Xo(w) -bxo(w) 


1 -b VcXi.^ = 0) 
Xo(w = 0) 


X(w), (6) 
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where xo is the independent-particle polarizability. We 
evaluate this for a real material with a long-range 
Coulomb interaction and ab initio band structure. 
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FIG. 2. (Color online) Optical absorption spectrum of LiF. 
Continuous line (red): physical solution; dashed line (blue): 
unphysical solution; Dots (black): experiment [iflil . The ver¬ 
tical line indicates the position of the quasiparticle gap. 


Again, this equation has two solutions that can be ob¬ 
tained by two iteration schemes, analogous to those in 
Eq. (IH) [^. Figure [2] shows the absorption spectrum 
of LiF obtained with the two iteration schemes, as well 
as the experimental result. The experimental spectrum 
shows a strongly bound exciton with a binding energy 
of approximately 1.4 eV 0. The spectrum obtained 
from iterating Eq. m with the analogue of scheme (I) is 
qualitatively correct, since it also shows a strongly bound 
exciton. The remaining discrepancies with respect to ex¬ 
periment are due to the approximate form of /xc . In¬ 
stead, iterating Eq. ([5]) within scheme (II) makes the ex¬ 
citon disappear completely. This means that, in absence 
of experimental results, one risks to run into a wrong 
solution, which would make the theory non-predictive. 
However, as we showed, the problem can be overcome 
using the appropriate iteration scheme (I). Note that op¬ 
tical properties can also be calculated from e = 1 — VcP, 
where the irreducible polarizability P obeys a Dyson 
equation with kernel fxc = 1/[(1 ~ VcP)xo]- In this case, 
the appropriate scheme is scheme (II). This is explained 
in the supplemental material Q. 

So far we have looked at the map Gq ^ G (and 
Xo x). We now focus on the inverse map Gq •<— G. 
This map is needed in problems of embedding, where one 
optimizes an auxiliary quantity Go in order to produce 
certain properties of a real system (contained in G). The 
inverse map is also crucial when one wants to express 
a functional in terms of dressed instead of bare quanti¬ 
ties. The most prominent example is the Luttinger-Ward 
(LW) functional, where the self-energy is given in terms 
of G instead of Go [2, 21-^. For the LW functional to 
be properly defined, the map Go •<— G should be unique. 

Within the 0PM, consider a system with the bare 
Green’s function yg, and with the exact, interacting 


Green’s function y given by Eq. (Ill)- We now fix y 
and examine whether the inverse map zo •<— ?/ unam¬ 
biguously leads to zo = Vo- With the exact self-energy 
s[zo] = — ^uzo of Eq. (jT]), the exact Dyson equation of 
this problem reads 

zo = y+ i^uyzl, (7) 

in which y is known and zq is to be determined. This 
equation has again two solutions: 




2 + V ±^{2-V)^ 
2V 


( 8 ) 

where Zq = zq/ yo and we used Eq. ([T|) . The square root 
in Eq. dS]) equals the absolute value |2—Ej. Because 2—V 
changes sign at E = 2, the physical solution Eo = 1 is 
obtained by Zq for E < 2 and by Zq for E > 2 (see 
Fig.p. In other words neither of the two solutions gives 
Zq = Yq for all E. As a consequence one has to change 
sign in front of the square root in Eq. (|5]) at E = 2. This 
has important consequences for the iterative solution of 
Eq. (I7|): because scheme (I) yields the square root with 
positive sign, to obtain the map Go ^ G we need two 
different iteration schemes: one for 0 < E < 2 and the 
other for E > 2. This is different from the map Go —t G, 
where one solution gives the physical solution for all E, 
and hence a single iteration scheme suffices. 

The need to change iteration scheme is a serious prob¬ 
lem. Indeed, Kozik et al. Q pointed out that differ¬ 
ent iteration schemes, applied to Hubbard and Anderson 
models, lead to different solutions which cross at a cer¬ 
tain interaction. Our 0PM results provide the missing 
explanation: keeping the labels (A) and (B ) of 0, the 
two iteration schemes correspond to 
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(rt+1) 




(A), 

(9) 

(B). 

(10) 


We report the results in Fig.[T] Scheme (A) converges to 
the physical solution for E <2 and to the nonphysical 
solution for E > 2. Instead, scheme (B) converges to the 
nonphysical solution for 2.IZ< E < 2 and to the physical 


solution for 2 < E <6 [2^. These results are strictly 


analogous to those obtained by Kozik et al. for Hubbard 
and Anderson models. They can be understood from the 
fact that scheme (A) creates a continued fraction with 
positive square root, whereas in scheme (B) the sign of 
the continued fraction is changed. 

This sign problem is a priori a disaster because there 
is no unique prescription of how to avoid unphysical so¬ 
lutions. The 0PM highlights the reducible polarizabil¬ 
ity 0 


X o 2-E 
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( 11 ) 
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FIG. 3. (Color online) One-point model (OPM): y as a func¬ 
tion of the interaction V for 4 <V < 100. Solid line (black): 
exact solution; dotted line (blue): Hartree-Fock (HF); dashed 
line (red): strong-interaction HF (SIN-HF). Inset: V as a 
function of V for 0 < H < 5. 


as critical quantity that changes sign at the crossing V = 
2 [^. At the same time, for y > 2 the perturbation 
expansion of y, in Eq. o, diverges. This is in line with 
Ref. 0 in which a breakdown of perturbation theory is 
linked to an eigenvalue of the polarizability that crosses 
zero, becoming negative. Our result confirms that one 
should inspect the reducible polarizability as a function of 
the interaction to detect problems of perturbation theory. 


Inverting a map between functionals requires a careful 
definition of their domain The multiple solutions 

are the price to pay for the fact that we have not con¬ 
sidered this definition in the above discussion. This can 
be understood as follows: if there were two solutions for 
Go ■«— G, one could obtain the same dressed Green’s func¬ 
tion from two different Gq and, hence, from two different 
external potentials. Since the diagonal of G is the density, 
the Hohenberg-Kohn theorem [23 states that there can 
only be one external potential, and hence one Go, cor¬ 
responding to each G. This means that any additional 
solution Go is unphysical, in the sense that it cannot be 
constructed from the solution of a one-body Schrodinger 
equation. Equivalently, it cannot be written as a sum of 
simple poles, each with a strength normalized to one. By 
imposing this condition, one can therefore eliminate un¬ 
physical solutions. In the OPM this trivially corresponds 
to the requirement Zq = 1, which alrea^ implies the so¬ 
lution. In the supplemental material we work out a 
more realistic case where the definition of the domain de¬ 
fines the solution. It should be noted that when Go is an 
embedding Green’s function the discussion is more com¬ 
plicated, because one searches for a fictitious Go with a 
frequency dependence that can differ from that of a Go 
resulting from a static potential. 

With the map Go ^ G one can construct the exact 


self-energy as a functional of G. Using Eq. ([5]) in the 
exact self-energy given in Eq. m we obtain 

s^[y,u\ = (l± v^l -2uy2^ (12) 

1 1 1 

= -Yy^Yy^r 

(13) 



The Dyson equation with the two self-energies of Eq. m 
leads to two different Green’s functions: the physical so¬ 
lution given in Eq. m is obtained from s for U < 2 
and from s’*' for U > 2, and y = 0 from s’*' for any 
V. Therefore, for weak interaction using the exact self¬ 
energy one obtains only one solution, the physical one, 
contrary to, e.5., the HF approximation. We note that at 
the point where and s~ meet (at U = 2) the deriva¬ 
tive ds^/dy diverges. This could explain the divergence 
of SY^/SG observed in Ref. [3] for a 2D Hubbard model. 
Note that this divergence occurs at the point where one 
of the eigenvalues of the polarizability crosses zero. 

In Eq. (fUTll we Taylor expanded the square root. The 
convergence radius is infinite, since 0 < 2uy^ < 1, as can 
be shown using Eq. ©• Interestingly, the sum of the 
first two terms in (1131) (upper sign) is the first term of 
an expansion of the self-energy for strong interaction . 
The remaining terms constitute an expansion in terms 
of a quantity that is proportional to u and converges 
for all physical y. This means that one can use pertur¬ 
bation theory over the whole interaction range , bu t in 
two different ways for the two different regimes [29j. To 
lowest order, this corresponds to HF, = —\uy, for 
weak interaction, and gSiN-HF _ _i_gHF^ strong in¬ 
teraction. We call this functional strong-interaction HF 
(SIN-HF). Both self-energies yield two solutions. We re¬ 
port the physical solution for these two approximations 
in Fig. 131 While HF clearly fails for strong interaction, 
SIN-HF is exact in the strong interaction limit and per¬ 
forms well for U > 4, while it is worse than HF for 
U < 4. It is important to note that the physical SIN- 
HF solution is obtained for U >1 with the iteration 
scheme I/y(”+i) = -|- — 1. Indeed, as 

it is also clear from the example of TDDFT, and dis¬ 
cussed in the supplemental material Q, the appropriate 
iteration scheme depends on the formulation of the prob¬ 
lem. We suggest the OPM as a powerful tool to examine 
which scheme one should use for a given problem and 
interaction range. 

In conclusion, we have demonstrated that with a sim¬ 
ple but general one-point model one can understand 
and solve structural problems of many-body perturbation 
theory. In particular, one can use it sort out the multi¬ 
ple solutions of the non-linear Dyson equation by choos¬ 
ing the appropriate iteration scheme. We have shown 
that for the map Go —>■ G a single iteration scheme suf¬ 
fices to obtain the physical solution for all interaction 
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strenghts. Instead, for the inverse map Gq G one has 
to change iteration scheme at the interaction strength 
at which the reducible polarizability changes sign and 
perturbation theory of G in terms of Gq starts to di¬ 
verge. Nonetheless, we have proved that even for strong 
interaction one can use a perturbative expression for the 
self-energy in terms of G, which differs from the usual 
LW functional. By presenting analogous results for real 
systems, and by comparing with numerical results in the 
literature, we have shown that these conclusions go far 
beyond the one-point model. We expect that they will 
have an impact on many other questions in the domain 
of many-body physics. 
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